using Pkg
Pkg.activate(".")

using DataFrames, CSV, StatsPlots, Downloads, DataFramesMeta, Dates, Statistics, Turing, Memoization, ReverseDiff, Serialization

Firearms and Violent Crime, recent trends

if !Base.Filesystem.ispath("data/nics-firearm-background-checks.csv")
    Downloads.download("https://github.com/BuzzFeedNews/nics-firearm-background-checks/raw/master/data/nics-firearm-background-checks.csv","./data/nics-firearm-background-checks.csv")
end

if ~Base.Filesystem.ispath("data/fipscodes.csv")
    Downloads.download("https://www2.census.gov/geo/docs/reference/state.txt","./data/fipscodes.csv")
end
fipscodes = @chain CSV.read("data/fipscodes.csv",DataFrame) begin  
    @subset(:STATE .< 60)
end


# manually downloaded CDC wonder gun homicide data using online form 
gunhomicidesnew = CSV.read("./data/wonder-gun-homicide-byyear.csv",DataFrame)
gunhomicidesold = CSV.read("./data/cdc-wonder-firearm-homicide-1979-1998.csv",DataFrame)
gunhomicides = [gunhomicidesnew ; gunhomicidesold]
gunhomicides = @chain rename(gunhomicides,Dict("Crude Rate" => "crudedeathrate", "State Code" => "StateCode", "Year" => "year", "Year Code" => "YearCode")) begin
    @subset(.! ismissing.(:Deaths))
    @transform(:crudedeathrate = map(x -> isnothing(x) ? missing : x, tryparse.(Float64,String.(:crudedeathrate))),
        :state = String31.(:State))
    @subset(.! ismissing.(:crudedeathrate))
    @orderby(:StateCode,:year)
end


gunchecks = CSV.read("./data/nics-firearm-background-checks.csv",DataFrame)
startpops = @chain gunhomicides, @subset(:year .== 1999), @select(:State, :StartingPop = :Population,:StartDeathRate = :crudedeathrate)



gunsalesest = @chain gunchecks,
    @orderby(:state,:month),
    @select(:state,:month,:salesest = 1.1 .* (:handgun .+ :long_gun) .+ 2.0 .* :multiple),
    groupby(:state),
    @transform(:cumsales = cumsum(:salesest),:year=year.(:month)),
    @subset(month.(:month) .== 6)


joined = leftjoin(gunhomicides,gunsalesest, on=[:State => :state,:year => :year]; matchmissing=:notequal,makeunique=true)

joined = @chain leftjoin(joined,startpops; on=[:state => :State],makeunique=true) begin
    @subset(.! ismissing.(:Deaths) .&& .! ismissing.(:StartingPop))
    @subset(:StateCode .< 60)
    @transform(:cumgunrate = (:cumsales .+ 1.0 .* :StartingPop) ./ :Population )
end

joined = @orderby(leftjoin(joined,fipscodes,on = :state => :STATE_NAME),:StateCode,:year)


display(plot(joined.cumgunrate,joined.crudedeathrate; group=joined.state, legend=false, xlab="Gun Rate (guns/capita)",
    ylab="Firearm Homicide Rate (deaths/100kcapita/yr)", title="Gun Homicide vs Gun Prevalence", alpha=.5))

display(plot(joined.year,joined.crudedeathrate; group=joined.state, legend=false, xlab="Year",
    ylab="Firearm Homicide Rate (deaths/100kcapita/yr)", title="Gun Homicide vs Time", alpha=.5))

How has "Constitutional Carry" affected violence in states that pass those laws? Let's take a first look at the data, by plotting the firearms homicide rate through time, with a marker for the onset of Constitutional Carry as found in the Wikipedia article on Constitutional Carry.

concarry = leftjoin(fipscodes,CSV.read("./data/ConstCarryDates.csv",DataFrame),on = :STUSAB => :StateCode)
concarry.date = map(x -> ismissing(x) ? missing : Date(div(x,10000),div(mod(x,10000),100),1),concarry.LawDate)
concarry.year = map(x -> ismissing(x) ? missing : year(x),concarry.date)
concarry = @subset(concarry,:STATE .<= 59) #ignore minor territories

let p = []
    for (st,stusab,statename,statens,lawdate,stdate,year) in eachrow(concarry)
        d = year
        sub = @subset(joined,:state .== statename)
        if(nrow(sub) > 0)
            pp = plot(sub.year,sub.crudedeathrate,ylim=(0,15),size=(500,500),
                title="$statename death rate vs time",xlab="Year",ylab="Gun Homicide/100k/yr",legend=false)
            if(!ismissing(d))
                pp = plot!([(d,0),(d,15)])
            end
        push!(p,pp)
        end
    end
    plot(p...; size = (2000,2000))
end

States by Geographic Group

Looking at the raw data we can see that there are a number of states that have an "accelerating" trend of increased firearms homicide violence in the period from about 2015 onward. Some of these have converted to Constitutional Carry near the beginning of that period or during that period, with the trend beginning before the carry provision for several of the states. Also many states have passed their laws in the last few years and have no data post-law.

Let's look at these states grouped into geographic regions, because trends in regions may be more obvious. We will start with the following groupings:

stategroups = [["WA","OR","CA","NV","AZ"],["ID","MT","WY","UT","CO","NM","TX"],["ND","SD","NE","KS","OK","IA","MO"],
    ["MN","WI","IL","IN","MI","OH"],["AR","LA","MS","AL","GA","FL","SC"],["TN","NC","KY","WV","VA","MD","DE"],
    ["PA","NJ","NY","CT","RI","MA","VT","NH","ME"],["AK","HI","DC"]]

stategroupdf = DataFrame(STUSAB = reduce(vcat,stategroups),stgroup = reduce(vcat,[[k for j in 1:length(stategroups[k])] for k in 1:length(stategroups)]))
stategroupdf = leftjoin(fipscodes,stategroupdf; on = :STUSAB)

@assert(length(unique(reduce(vcat,stategroups))) == 51) # 50 states plus DC

for statelist in stategroups
    sub = joined[in.(joined.STUSAB, Ref(statelist)),:]
    ccdates = @chain @subset(concarry, in.(concarry.STUSAB,Ref(statelist))) begin
        @subset(.! ismissing.(:year))
    end
    
    ccdatesenact = leftjoin(ccdates,sub, on = [:STUSAB => :STUSAB, :year => :year], makeunique = true)
    ccdatesenact = @subset(ccdatesenact, .! ismissing.(:crudedeathrate))
    #display(ccdatesenact)
    #display(sub)
    p = plot(sub.year,sub.crudedeathrate,group=sub.STUSAB,ylim=(0,20),legend=:topleft,linewidth = 3, thickness_scaling=1,
        title="Firearm Homicide Rate vs Time", ylab="Homicides/100k population/yr")
    if nrow(ccdatesenact) > 0
        p = scatter!(ccdatesenact.year,ccdatesenact.crudedeathrate,markersize=6,group = ccdatesenact.STUSAB)
    end
    display(p)
end

States by Relative Change

Each state has its own overall level of firearm violence, but often trends in the same direction relative to its region even if its overall level is higher or lower. This suggests a model in which we measure each state using a dimensionless number relative to some baseline level. For example the average rate in the years 1999,2000,2001

baselines = @chain joined begin
    @subset(in.(:year,Ref((1999,2000,2001))))
    @by(:STUSAB, :baseline = mean(:crudedeathrate))
end

bljoined = @chain leftjoin(stategroupdf,leftjoin(baselines,joined,on = :STUSAB, makeunique=true),on=:STUSAB,makeunique=true) begin
    @subset(.! ismissing.(:crudedeathrate))
    @orderby(:stgroup,:STUSAB,:year)
end

function getccd(gr)
    ccd = innerjoin(gr,@subset(concarry,.! ismissing.(:year)), on = [:STUSAB => :STUSAB, :year => :year], makeunique = true)
    ccd = @subset(ccd, .! ismissing.(:crudedeathrate))
    ccd
end

function makerelplot(gr)
    ccd = getccd(gr)
#    display(ccd)
    p = plot(gr.year,gr.crudedeathrate ./ gr.baseline,group=gr.STUSAB,legend=:topleft,ylim=(0.0,2.5),xlab="Year",ylab="Relative Rate") 
    if nrow(ccd) == 1
        scatter!(ccd.year,ccd.crudedeathrate ./ ccd.baseline)
    elseif nrow(ccd) > 1 
        scatter!(ccd.year,ccd.crudedeathrate ./ ccd.baseline,group = ccd.STUSAB)
    end
    p
end

makerelplot (generic function with 1 method)

What has been the relative change?

Each state plotted in a group with other geographically similar states. Firearm homicide rate relative to the recorded rate in the year 1999.

plots = [makerelplot(gr)
    for gr in groupby(bljoined,:stgroup)]
    
plot(plots... ; size=(1000,1000),linewidth=2,xlab="Year",ylab="Relative Rate")

Modeling the process with Bayesian Models

We can try to estimate the effects of these laws by building a model of the overall process. Because we can't resurrect people from the dead, the homicide rate can never be negative. It makes sense then to model the homicide rate on a logarithmic scale.

We estimate the behavior of the states as an overall level, plus a shape that is common to the region. Individual states are then allowed a small perturbation to the regional shape.

The effect of constitutional carry law would then be estimated as the log firearm homicide rate which exceeds the state level counterfactual estimate in the years after the passage of the law.

We impose an informal smoothness requirement on the counterfactual estimates by using a compact radial basis function expansion with one center every 4 years, and a maximum radius of 7 years so that there is overlap between adjacent centers. 

function bumpfun(x,c,scale)
    stdx = (x-c)/scale
    if stdx < -1.0 || stdx > 1.0
        0.0
    else
        exp(1.0-1.0/(1.0-stdx^2)) # goes to zero at -1 and 1, and 1 at x=0
    end
end

function timeser(x,coefs,centers,scale)
    f = 0.0
    for (a,c) in Iterators.zip(coefs,centers)
        f += a*bumpfun(x,c,scale)
    end
    return f
end

function laweffect(yr,rate,start)
    if yr >= start
        1.0-exp(-rate*(yr-start))
    else   
        0.0
    end
end


statecoefpert = log(1.05)/2.0

include("model.jl") ## load the model, this is a separate file so we can save the output unless it changes

modeljoined = @chain leftjoin(joined,stategroupdf; on = :STUSAB, makeunique=true) begin
    @subset(.! ismissing.(:crudedeathrate))
end

lawdates = @chain leftjoin(DataFrame(STATE=1:55),concarry; on=:STATE) begin
    @subset(:STATE .< 60)
    @rtransform(:year = ismissing(:year) ? 3000 : :year)
    @orderby(:STATE)
end

centers = [i for i in (minimum(modeljoined.year)-4):4:(maximum(modeljoined.year)+4)]
width = 7.0


modl = guns(modeljoined.StateCode,maximum(modeljoined.StateCode),modeljoined.stgroup,maximum(modeljoined.stgroup),
        centers, width, modeljoined.year,log.(modeljoined.crudedeathrate),lawdates.year,statecoefpert)

setadbackend(:reversediff)
Turing.setrdcache(true)

savedfile = "./saved/samples.dat"
global s = []
if stat(savedfile).mtime > stat("./model.jl").mtime
    global s = deserialize(savedfile)
else
    global s = sample(modl,NUTS(500,0.8),MCMCThreads(),500,3)
    serialize(savedfile,s)
end

Having sampled the model, we can then compute summary graphs of the results. We will show the posterior density for the coefficient of the magnitude of the law effect for each state separately.

Also we will plot the actual log(homicide rate) together with the state specific counterfactual estimates.

lawco = group(s,"statelawcoef")

statenames = Dict()
for (st,stusab) in Iterators.zip(lawdates.STATE,lawdates.STUSAB)
    statenames[st] = stusab
    end

function plotlawts(s,lawdates,modeljoined)
let pl = []
    lawco = group(s,"statelawcoef")

    for i in @select(@subset(lawdates,:year .< 2022),:STATE).STATE
        den = density(s[:,Symbol("meanlawcoef"),:] .+ lawco[:,Symbol("statelawcoef[$i]"),:], title = "State $i = $(statenames[i])",xlim=(-log(2),log(2)))
        plot!([(0.0,0.0),(0.0,5.0)], color="red",ylim=(0.0,5.0))
        push!(pl,den)
    end
    println("State law coefficient values")
    display(plot(pl...; size = (1000,1000)))

    display(density(s[:,:lawrate,:], title= "Law effect onset rate (1/yr)"))

    modeljoined.logdrate = log.(modeljoined.crudedeathrate)

    regions = Dict()
    for r in eachrow(stategroupdf)
        regions[r.STATE] = (group=r.stgroup,code=r.STUSAB)
    end
    pl = []
    samps = sample(1:500,10)

    for st in unique(modeljoined.STATE)
        sub = @subset(modeljoined, :STATE .== st)
        region = regions[st].group
        p = plot(sub.year,sub.logdrate,linewidth=3,title="$(st) = $(regions[st].code)",ylim=(-0.5,2.75))
        push!(pl,p)
        for samp in samps
            statelawcoef = s[samp,Symbol("statelawcoef[$st]"),1]
            statecoefs = [s[samp,Symbol("statecoefs[$st][$i]"),1] for i in 1:length(centers)]
            statebase = s[samp,Symbol("statebase[$st]"),1]
            regioncoefs = [s[samp,Symbol("regioncoefs[$region][$i]"),1] for i in 1:length(centers)]        
            
            lawrate = s[samp,:lawrate,1]
            startlaw = lawdates[lawdates.STATE .== st,:].year[1]
            startlaw = ismissing(startlaw) ? 3000 : startlaw
            years = 1979:2020
            pred = [statebase + 
                timeser(yr,regioncoefs .+ (statecoefpert .* statecoefs),centers,width) for yr in years]
            plot!(years,pred; color="orange",alpha=0.5)
            if startlaw < 3000
                plot!([(startlaw,-10.0),(startlaw,10.0)],color="red")
            end
            regpred = [statebase + 
                timeser(yr,regioncoefs ,centers,width) for yr in years]

        end
    end
    println("Each state shown with its locally estimated counterfactual (orange), actual (dark blue), and date of CC law if any (red vertical)")
    display(plot(pl...; size =(2000,3000), legend=false))
end

density(s[:,:meanlawcoef,:],title="Mean law coefficient")
end

plotlawts(s,lawdates,modeljoined)

State law coefficient values
Each state shown with its locally estimated counterfactual (orange), actual
 (dark blue), and date of CC law if any (red vertical)

Shall-issue vs permitless

One plausible explanation for the relative lack of clear measurable effect from permitless carry laws is that permitless actually may have had very little effect on how many people carry firearms in public. Many of these states were shall-issue and many people who wanted to carry may have had CCW permits issued before the permitless law passed. Therefore the change in number of people carrying may have been minimal.

Using the date on which states converted to shall-issue permits, with a similar analysis, we can determine whether more of an effect was visible after that legal change. Some states converted from no-issue to shall issue, while others converted from may-issue to shall-issue.

# dates collected from https://en.wikipedia.org/wiki/History_of_concealed_carry_in_the_United_States
shallissdates = CSV.read("data/ShallIssueDates.csv",DataFrame)

shallissue = leftjoin(fipscodes,shallissdates; on = :STUSAB)

silawdates = @chain leftjoin(DataFrame(STATE=1:55),shallissue; on=:STATE) begin
    @subset(:STATE .< 60)
    @rtransform(:year = ismissing(:ShallIssueDate) ? 3000 : :ShallIssueDate)
    @orderby(:STATE)
end

modl2 = guns(modeljoined.StateCode,maximum(modeljoined.StateCode),modeljoined.stgroup,maximum(modeljoined.stgroup),
        centers, width, modeljoined.year,log.(modeljoined.crudedeathrate),silawdates.year,statecoefpert)



savedfile = "./saved/shallisssamples.dat"
global sisam = []
if stat(savedfile).mtime > stat("./model.jl").mtime
    global sisam = deserialize(savedfile)
else
    global sisam = sample(modl2,NUTS(500,0.8),MCMCThreads(),500,3)
    serialize(savedfile,sisam)
end

Shall Issue coefficients

plotlawts(sisam, silawdates,modeljoined)

State law coefficient values
Each state shown with its locally estimated counterfactual (orange), actual
 (dark blue), and date of CC law if any (red vertical)

The Take Home

Many states that have passed constitutional carry laws have done so recently enough that there is no data on homicides available from CDC Wonder yet.

For those which passed the law long enough ago to have some post-law data, still many of them were in the last few years leaving only a few years of post-law data. The last few years is a time period where homicides have increased nearly nationwide and it is challenging to estimate the effect of a law when it is imposed on top of a nonlinearly changing trend. Estimates of a law effect created by comparing the actual data to counterfactuals estimated as small perturbations to the regional trend lead to estimates of the average law coefficient that has uncertain sign, with a most likely value of about 0.05 or so, meaning that gun homicides may increase by a factor of around 1.05, however the ability to resolve this number is so poor that the number could be anything from about -0.1 to +0.25 (multiplying the crime rate by anywhere from 0.9 to 1.28.

Evidence for a consistent average net positive or negative effect of passing constitutional carry laws simply doesn't exist in the data so far.

In states such as AK, AZ, and AR which passed their laws long enough ago that we might expect to be able to estimate the effect, the average effect appears to be close to 0.0 with uncertainty relatively symmetrically distributed to both sides. The strongest estimates we have are near zero.

KY, and MS appear to be outliers, but even in these cases it's ambiguous whether the law resulted in an increase or a decrease, even though increase appears to be the more probable direction.

The evidence simply does not show a consistent effect increasing or decreasing firearm homicides after liberalization of concealed carry laws.

Published from analyze.jmd using Weave.jl v0.10.10 on 2022-07-08.